3 Algoritmo
set.seed(288) ## setando a semente para gerar sempre os mesmos resultados da amostra
amostra <- rgamma(112,420,10) # gerando amostra da gammaA amostra segue uma distribuição \(gama(\alpha = 420,\beta = 10)\), logo esses são os verdadeiros parâmetros que tentaremos estimar através do método bayesiano, note que \(E(Y) = \frac{\alpha}{\beta} = 42\), logo nesse exemplo criado a empresa em média não cumpre com o plano contratado de 50 Megabits. Obs: Para a análise vamos considerar todos os parâmetros desconhecidos e seus valores reais não serão mencionados a partir daqui.
n = length(amostra) # tamanho da amostra
media = mean(amostra) # media amostral
variancia = var(amostra) # variância amostral
soma_y = sum(amostra) # soma da amostra
soma_log_y <- sum(log(amostra)) # soma do log da amostra
media ## valor da media## [1] 42.09928Aqui foi calculado algumas estatísticas que serão utéis para mais adiante, é possivel notar que a média foi bem consistente com seu verdadeiro parâmetro, indicando uma boa qualidade da amostra.
3.1 Escolha das Prioris
Como o objetivo é averiguar se a operadora está ofertando uma conexão conforme o previsto, se escolheu prioris para alpha e beta de tal maneira que a operadora não possa atribuir o resultado negatativo do teste às prioris. Abaixo foi mostrado esse Gráfico de uma \(gama(\alpha = 50,\beta = 1)\) e perguntado se representa bem a variabilidade da internet, João concordou com a distribuição. A distribuição segue a premissa de que a operadora distribui em média uma velocidade de conexão mínima, portanto com esperança igual a 50, logo a operadora não pode reclamar de viés da priori.
alpha <- 50
beta <- 1
x <- seq(qgamma(0.0001,shape = alpha,rate = beta),qgamma(0.9999,shape = alpha,rate = beta),0.01)
p <- plot_ly(x = x, y = dgamma(x,shape = alpha, rate = beta), type = 'scatter', mode = 'lines', fill = 'tozeroy')
p## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.Logos as prioris \(\alpha_0\) e \(\beta_0\) foram escolhidas de tal maneira que \(E(\alpha_0) = 10\) e \(E(\beta_0) = 1/5\) e com variâncias bem altas.
# constante para alpha0, note que quando menor k1 maior será a variancia!
k1=1/5
## Priori para alpha
lambda_alpha = 50*k1 # lambda da priori alpha
nu_alpha = k1 # nu da priori alpha
x <- seq(0,qgamma(0.995,shape = lambda_alpha,rate = nu_alpha),0.01)
p <- plot_ly(x = x, y = dgamma(x,shape = lambda_alpha, rate = nu_alpha), type = 'scatter', mode = 'lines', fill = 'tozeroy')
p## Priori para beta
k2 = 1/2
lambda_beta = k2 # lambda da priori beta
nu_beta = k2 # nu da priori betax <- seq(0,qgamma(0.99,shape = lambda_beta,rate = nu_beta),0.0001)
p <- plot_ly(x = x, y = dgamma(x,shape = lambda_beta, rate = nu_beta), type = 'scatter', mode = 'lines', fill = 'tozeroy')
pPelos gráficos conseguimos perceber que as prioris são pouco informativas como queriamos.
### Passo 1
n_simu <- 300000 # número de simulações
simu <- matrix(nrow = n_simu,ncol = 3) # inicializando matriz com as simulações
colnames(simu) <- c("alpha","beta","x") # nome das colunas da matriz
#simu[1,1:3] <- c(media^2/variancia,media/variancia,media) # calculando a primeira linha
simu[1,1:3] <- c(1000,100,media) # calculando a primeira linhavar2 = 10 ## variancia da gama que servira como distribuicao proposta
simular_proposta <- function(mu){ ## gera valor de uma gamma com media igual a mu e variancia igual a var2
rgamma(1,shape = mu^2/var2,rate = mu/var2)
}
proposta <- function(y,mu){ ## calcula a densidade da distribuicao proposta no ponto y dado mu
dgamma(y,shape = mu^2/var2,rate = mu/var2)
}priori_alpha <- function(y){ ## calcula a densidade para a priori de alpha
dgamma(y,lambda_alpha,nu_alpha) + (dgamma(y,lambda_alpha,nu_alpha) == 0)*.000001 ## esse ultimo termo é só para valores de priori igual a 0 não dar problema com log
}
log_vero <- function(alpha,beta){ ## log verossimilhança
n*alpha*log(beta) - n*lgamma(alpha) + (alpha-1)*soma_log_y
}## só para calcular o indice de aceitação de alpha
aceita = c()
set.seed(2465)
for (i in 2:n_simu) {
## pegando os valores anteriores da cadeia
alpha_anterior <- simu[i-1,1]
beta_anterior <- simu[i-1,2]
### Metropolis Hasting para gerar alpha
alpha_proposto <- simular_proposta(mu = alpha_anterior) ## simulando um valor proposto para alpha da distribuicao gama com media igual ao alpha do passo anterior
## calculo do da log_verossimilhanca no ponto (alpha proposto,beta do passo anterior)
loglik.atual <- log_vero(alpha_proposto,beta_anterior)
## calculo do da log_verossimilhanca no ponto (alpha passo anterior,beta do passo anterior)
loglik.anterior <- log_vero(alpha_anterior,beta_anterior)
## calculo da densidade da log priori no ponto alpha proposto/alpha anterior
log.prior.alpha.atual <- log(priori_alpha(alpha_proposto))
log.prior.alpha.anterior <- log(priori_alpha(alpha_anterior))
## calculo da densidade da log da distribuição proposta no ponto do alpha proposto dado o alpha anterior e do alpha anterior dado o alpha proposto
log.proposta.atual <- log(proposta(alpha_proposto,alpha_anterior))
log.proposta.anterior <- log(proposta(alpha_anterior,alpha_proposto))
logpi.atual <- loglik.atual + log.prior.alpha.atual - log.proposta.atual
logpi.anterior <- loglik.anterior + log.prior.alpha.anterior - log.proposta.anterior
u <- runif(1)
aceita[i-1] <- log(u) <= (logpi.atual - logpi.anterior)
alpha_atual <- aceita[i-1]*alpha_proposto + (1- aceita[i-1])*alpha_anterior
###############
# Atualizar Beta
###############
## gerando o valor de beta pela posteriori condicional
beta_atual <- rgamma(1, alpha_atual*n + lambda_beta, rate = soma_y + nu_beta)
## gerando um valor de gamma (desnecessário, mas quis deixar :P)
x_atual <- rgamma(1, alpha_atual,beta_atual)
## Guardando os valores
simu[i,] <- c(alpha_atual, beta_atual,x_atual)
}Taxa de aceitação da cadeia
## [1] 0.45337153.2 Resultado da cadeia
Com esse primeiro gráfico podemos notar que a cadeia convergiu, apesar que tenha demorado um pouco.
Retirando as 5000 primeiras observações da cadeia podemos perceber que não é mais necessário retirar mais observações
No gráfico de autocorrelação podemos perceber um grande, as observações são muito correlacionados e só começam a diminuir com um número muito grande de espaçamento entre as observações, com aproximadamente 1500 observações de espaçamento as observações se tornam independentes, o problema é que com isso precisamos de 1500 simulações para gerar um valor da posteriori independente, que causa um enormo custo computacional.
ps: Nos tentamos de muitas maneiras de diminuir esse problema, mas nenhuma foi efetiva :(
3.3 Interpretação dos Resultados
##
## Iterations = 5000:299000
## Thinning interval = 1500
## Number of chains = 1
## Sample size per chain = 197
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## alpha 210.735 25.2940 1.80212 1.80212
## beta 5.004 0.6041 0.04304 0.04304
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## alpha 163.569 192.89 210.279 228.851 258.058
## beta 3.868 4.59 4.983 5.429 6.123A esperança para \(\alpha\) tem valor estimado igual a 210.735 e mediana igual a 210.279, o intervalo de credibilidade central de 0.95 de probabilidade é dado por (163.569;258.058). Para \(\beta\) a esperança é igual a 5.004 e mediana igual a 4.983 com intervalo de credibilidade central de 0.95 de probabilidade dado por (3.868;6.123).
Se considerarmos a esperanças dos parâmetros então \(E(Y)= \frac{\alpha}{\beta} = 42.11331\), como esse valor envolve dois parametros é dificil gerar um valor para essa estimativa, talvez fosse o caso de simular com uma noma parametrização.
3.4 Gerar com outra distribuição proposta para \(\alpha\)
### Passo 1
n_simu <- 300000 # número de simulações
simu2 <- matrix(nrow = n_simu,ncol = 3) # inicializando matriz com as simulações
colnames(simu2) <- c("alpha","beta","x") # nome das colunas da matriz
#simu[1,1:3] <- c(media^2/variancia,media/variancia,media) # calculando a primeira linha
simu2[1,1:3] <- c(1000,100,media) # calculando a primeira linhasimular_proposta2 <- function(mu){ ## gera valor de uma qui quadrado com media igual a mu
rchisq(1,df=mu)
}
proposta2 <- function(y,mu){ ## calcula a densidade da distribuicao proposta no ponto y dado mu
dchisq(y,df = mu)
}## só para calcular o indice de aceitação de alpha
aceita = c()
set.seed(2465)
for (i in 2:n_simu) {
## pegando os valores anteriores da cadeia
alpha_anterior <- simu2[i-1,1]
beta_anterior <- simu2[i-1,2]
### Metropolis Hasting para gerar alpha
alpha_proposto <- simular_proposta(mu = alpha_anterior) ## simu2lando um valor proposto para alpha da distribuicao gama com media igual ao alpha do passo anterior
## calculo do da log_verossimilhanca no ponto (alpha proposto,beta do passo anterior)
loglik.atual <- log_vero(alpha_proposto,beta_anterior)
## calculo do da log_verossimilhanca no ponto (alpha passo anterior,beta do passo anterior)
loglik.anterior <- log_vero(alpha_anterior,beta_anterior)
## calculo da densidade da log priori no ponto alpha proposto/alpha anterior
log.prior.alpha.atual <- log(priori_alpha(alpha_proposto))
log.prior.alpha.anterior <- log(priori_alpha(alpha_anterior))
## calculo da densidade da log da distribuição proposta no ponto do alpha proposto dado o alpha anterior e do alpha anterior dado o alpha proposto
log.proposta.atual <- log(proposta(alpha_proposto,alpha_anterior))
log.proposta.anterior <- log(proposta(alpha_anterior,alpha_proposto))
logpi.atual <- loglik.atual + log.prior.alpha.atual - log.proposta.atual
logpi.anterior <- loglik.anterior + log.prior.alpha.anterior - log.proposta.anterior
u <- runif(1)
aceita[i-1] <- log(u) <= (logpi.atual - logpi.anterior)
alpha_atual <- aceita[i-1]*alpha_proposto + (1- aceita[i-1])*alpha_anterior
###############
# Atualizar Beta
###############
## gerando o valor de beta pela posteriori condicional
beta_atual <- rgamma(1, alpha_atual*n + lambda_beta, rate = soma_y + nu_beta)
## gerando um valor de gamma (desnecessário, mas quis deixar :P)
x_atual <- rgamma(1, alpha_atual,beta_atual)
## Guardando os valores
simu2[i,] <- c(alpha_atual, beta_atual,x_atual)
}
##
## Iterations = 5000:299000
## Thinning interval = 3500
## Number of chains = 1
## Sample size per chain = 85
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## alpha 213.54 24.4106 2.64770 2.64770
## beta 5.07 0.5869 0.06365 0.06365
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## alpha 172.006 194.555 210.408 228.851 261.990
## beta 4.085 4.607 4.999 5.446 6.229Nesse o gráfico de autocorrelacao ficou pior ainda, precisando de um espaçamento de 3500 observações.